In contrast from ideal solid, fluid, or gas, granular material has a unique behavior and can be treated as a model for studying statistics mechanics or condensed matter physics.^[1–3] Granular flow from hoppers has been widely studied during last few decades.^[4–9] Granular materials discharge constantly from hopper under gravity, which is different from normal fluids^[10,11] and the flow rate is independent of the hopper diameter D with some criteria satisfied.^[6] However, whether this constancy is related with Janssen effect is controversial.^[6,11,12] To quantitatively describe the flow rate, a widely accepted scaling formula called Beverlooʼs law was proposed by Beverloo et al. in 1961.^[4] This formula has a unified form for various types of hoppers: , where C is a dimensionless coefficient which depends on the material properties,^[11,13]R is the equivalent hydrodynamic radius, ρ_i is the equivalent density of the initial grain pack in the hopper, g is the acceleration of gravity, and A is the equivalent hydrodynamic section. This is valid when R is larger enough than the grains size.^[7,14] Recently, the experimental and numerical studies showed that this law is not exactly valid if the outlet size is very large and several modified formats of the formula were suggested.^[9,15,16]

The free fall arch (FFA),^[8,9,17] which is assumed to exist above the orifice, is useful for understanding the scaling law. In this assumption, the velocities of grains above the arch are negligible and below the arch grains fall freely under gravity.^[18] If the height of the FFA is proportional to the orifice diameter, then the grains velocities at the orifice should be proportional to .^[7] With A being the section that the grains flow through, the volume flowed out per unit time should be and the flow rate should be proportional to , where ρ_f is the equivalent density. However, recent results have shown that the assumption of FFA is not exactly valid as traditionally understood.^[18–20]

To decouple the relationship between the diameter of hopper D and the flow rate, two criteria were proposed formerly: one is and the other is .^[6] Here we focus on these criteria through large-scale numerical simulations. Discrete element method (DEM) is employed to simulate millions of spheres discharging from the hopper. In these simulations, the diameters of hoppers and orifices are changed in a large range to study the relationship between them and flow rates. Then, the deviation from Beverlooʼs law when the criteria are not fulfilled is investigated.

A considerable amount of work has studied granular flows using computer simulations^[8,21] and several approaches have been applied. To simulate millions of grains, this work was carried out on multiple GPUs by the DEM code we developed.^[22,23] In the DEM model, the interactions between the spheres are given by Hertz–Mindlin contact model.^[24,25] Supposing the spheres are identical, the normal and tangential contact forces between two contacted spheres are

Then, in the gravity field, the equations of motion of the spheres are

In this work, monosized glass spheres are randomly generated with a small volume fraction within three-dimensional (3D) flat bottomed hoppers, the material of which is the same as the spheres’. The number of spheres (30000–4360000) varies with different hopper sizes and thus packing heights could be guaranteed to be equal (about 330d, d is the sphere diameter). For each hopper diameter D, the diameter of rounded orifices D₀ varies from 6d to D − d. The spheres are inserted in the hopper space randomly with an initial volume fraction of 0.2 and then they are packed under gravity till the total kinetic energy of the packing is small enough ( ). Then, the orifices of the hoppers are opened to let the spheres flow. In view of the fluctuations during the flows, the flow rate in our simulations is calculated from the slopes in discharge profiles^[11] and cross-checked with parallel simulations using different initial packing constructions. The statistical values of velocities are obtained after being averaged over time. We further simulate the cases with different parameters. Table 1 shows the parameters of glass as a base case of our simulations.

Table 1.

Table 1.

Table 1. Parameters of the base case in our simulations. . Quantity Symbol Value Elastic modulus/GPa Y 72 Poissonʼs ratio ν 0.25 Friction coefficients μ 0.5 Diameter of spheres/m d 0.001 Density of spheres/ ρ 2500 Coefficient of restitution ε 0.4

Table 1. Parameters of the base case in our simulations. .

Figure 1 shows the flow rates of different hopper orifice diameters with hopper diameter D = 90d. The flow rates of different random initial packings are also plotted in this figure. From this figure, it can be seen that the flow rates become stable within 0.1 s even for the case that the orifice diameter is very close to the hopper size. There are a few fluctuations in the temporal profile of instantaneous flow rate, which is consistent with the previous study.^[27] This result also shows that the initial packing state does not influence the flow rate.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Temporal profiles of instantaneous flow rate with μ=0.5 and D = 90d. For parallel numerical experiments, 3 different initial packings were simulated.

The influence of the changes in material and geometrical parameters is studied by running a series of simulations. For each case, only one parameter differs from that of the base case while the others remain the same. The base case was introduced in the previous section. Figure 2 shows the flow rates of discharging hoppers with D = 120d and different parameters. When the orifice diameter is not large, all of the flow rate profiles present a unique scaling law and are well fitted by the Beverlooʼs law

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Flow rates plotted against D0 with different parameters when D = 120d. Only one parameter varies for each case while the other parameters are the same with those of the base case. The stars denote the flow rate when sphere diameter d=0.01 m.

Table 2.

Table 2.

Table 2. Fitted parameters for different cases. . Cases C k Base case 0.6 1.4 Y=7.2 GPa 0.61 1.4 Y=0.72 GPa 0.62 1.4 μ=0.2 0.68 1.4 μ=0.8 0.59 1.4 ε=0.8 0.61 1.4 ρ = 250 kg/m3 0.61 1.4 ρ = 25000 kg/m 3 0.61 1.4 d=0.01 m 0.61 1.4

Table 2. Fitted parameters for different cases. .

In following, large-scale hopper flows with various hopper and orifice diameters are simulated and three questions are raised and studied: 1) Is there a uniform scaling law for different hopper and orifice diameters? 2) What is the influence of D on the flow rate? and 3) Are the criteria good for Beverlooʼs scaling law?

Figure 3(a) plots the variations of flow rate of hopper with different hopper and orifice diameters. For each case of a certain hopper diameter, the flow rates are in good agreement with the 2.5 power law when the orifice diameter D is not very large, which indicates that there is a uniform scaling law for different hopper and orifice diameters. However, the discrepancy between flow rates of simulations and predicted values from Eq. (3) becomes larger as the orifice diameter D₀ approaches to the hopper diameter D (see Fig. 3(b)). Figure 4 shows the influence of D on the flow rate. As expected, when , as D increases, the flow rate decreases to a certain value, which shows the decoupling of hopper size and flow rate. For large D₀, the influence is obvious even if D=90d.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The relationship of flow rate with D0 for different hopper diameters with μ=0.5. (a) A log–log plot of flow rate; (b) a plot of flow rate ratio (the ratio of the simulation result and the predicted result using Eq. (3)) when D0 approaches to D.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Variations of flow rate ratio with increasing D.

To investigate the criteria for Beverlooʼs scaling law, two criteria that are introduced above are checked here and the results are shown in Fig. 5. For criterion 1, it is found that when , the flow rates can be well described by Eq. (3). When , there is a little discrepancy between the flow rates and prediction. For criterion 2, it is found that the flow rates can be well described by Eq. (3) when in our results.^[29] Interestingly, it indicates that there is still a scaling relation between the flow rate and orifice diameter if the ratio of orifice diameter and hopper size is fixed and less than 2 (see the line in Fig. 5(b)). Then the parameter C should be larger than that in Eq. (3). From the previous discussion, we can conclude that criterion 1 should be suggested as with large orifice diameter and criterion 2 should be suggested as . The new criterion 1 here is more stringent than that in Ref. [6], but the criterion 2 is looser.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Flow rate ratio plotted against with flow rate ratio varied with increasing D0 and (a) constant (b) constant .

Generally, the flow rates in hopper flows depend on the vertical velocities and volume fractions at the orifices. Figure 6 shows the similar spatial profiles of volume fractions in hopper flows with different D₀. It is found that in most space, the volume fraction is about 0.6. Above the orifices, the volume fraction is about 0.55 and the volume fraction decreases to 0.45 at the rim of the orifice. Meanwhile, from the ideal FFA assumption, the grains fall freely from a semi-spherical surface;^[8] but recently the free fall surface was suggested as a parabolic surface.^[9,30] Then, the vertical velocity at the orifice should be given by

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Spatial profile of volume fraction on vertical cross section with D = 60d and (a) D0 = 20d; (b) D0 = 30d; (c) D0 = 40d; (d) D0 = 50d.

The profiles of vertical velocity at the orifices are shown in Fig. 7 and are fitted by Eq. (4). From this figure, the results of D₀ = 6d, 10d, 20d, and 30d can be fitted to the equation, while the velocities at the large orifices (D₀ = 50d and 55d) fail. For the large orifice cases, the velocities in the central region are nearly the same and drop sharply at the orifice edges, which is similar to the profile of plug flow region.^[31] Figure 8 plots the spatial profiles of vertical velocity of the sections across the hopper symmetry axis. Within D₀ = 20d hopper, most of vertical velocities are very small and the vertical velocity increases rapidly above the orifice, which agrees with the FFA assumption.^[18,30] However, for the larger orifice, the grains flow rapidly within the hoppers and are accelerated a bit at the orifice, with only a very small stagnant zone at the bottom near the wall.^[32,33] The vertical velocity of grains above the free fall region cannot be neglected and the FFA assumption is invalid here. When D₀ is small, the variation of flow rate obeys Beverlooʼs law and when D₀ approaches D, the variation deviates from the law because of the change of the velocity field.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Spatial profile of vertical velocity on the central axis of the orifice. The lines are fitted results by Eq. (4) and the fitted values of vc are 0.22 m/s, 0.29 m/s, 0.46 m/s, and 0.55 m/s for D0 = 6d, 10d, 20d, and 30d, respectively.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Spatial profile of vertical velocity on vertical cross section with D = 60d and (a) D0 = 20d; (b) D0 = 30d; (c) D0 = 40d; (d) D0 = 50d.

In this paper, the granular flows discharged from a hopper with a large range of diameters of hopper and orifice are simulated and the criteria for decoupling of hopper size and flow rate are studied. When the hopper size is fixed, it is found that the Beverlooʼs law is valid except the orifice diameter is close to the hopper size. Then, the vertical velocities of grains above the free fall region are much larger and disobey the FFA assumption. The change of velocity field needs to be investigated quantitatively in the future. In this study, monosized spheres are simulated and the effect of dispersity should be considered for the next step.